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Abstract 

The present paper describes the results obtained from a one-dimensional time dependent 
numerical technique that simulates early flame propagation in a moderate to intense turbulent 
environment. Attention is focused on the development of a spark-ignited, premixed, lean 
methane/air mixture with the unsteady spherical flame propagating in homogeneous and 
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isotropic turbulence. A Monte-Carlo particle tracking method, based upon the method of 
fractional steps, is utilized to simulate the phenomena represented by a probability density 
function (PDF) transport equation. Gaussian distributions of fluctuating velocity and fuel 
concentration are prescribed. Attention is focused on three primary parameters that influence the 
initial flame kernel growth: the detailed ignition system characteristics, the mixture composition, 
and the nature of the flow field. The computational results of moderate and intense isotropic 
turbulence suggests that flames within the distributed reaction zone are not as vulnerable, as 
traditionally believed, to the adverse effects of increased turbulence intensity. It is also shown 
that the magnitude of the flame front thickness significantly impacts the turbulent consumption 
flame speed. Flame conditions studied have fuel equivalence ratio’s in the range (j) = 0.6 to 0.9 at 
standard temperature and pressure. 


Nomenclature 


Symbol 

Definition 

B 

pre-exponential factor 

B 

reduced pre-exponential factor 

C 

concentration 

c p 

Y p 

normalized mean fuel concentration, 

f7ct 

c'f 

r SI 

root mean square (rms) fluctuation of normalized fuel concentration 

c, 

l T 

coefficient of the integral length scale of turbulence 

D/Dt 

total or substantial derivative of 

Da 

Damkohler number 

Da 

turbulent Damkohler number 

E a 

activation energy (Joule/mole) 

E. 

ign 

actual energy deposited into gas by ignition source spark energy 

E 

ign-min 

minimum ignition energy 

E sp 

electrical energy in the capacitance spark 

k 

turbulent kinetic energy 

K 

flame stretch 

Ka 

Karlovitz number 
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Karlovitz flame stretch factor 

Karlovitz number for unstretched laminar flame conditions 

integral length scale of turbulence 

distance between spark electrodes 

Markstein number 

total number of particles in a cell 

number of particles per cell 

pressure, Joint Probability density distribution function 
radial coordinate 

( J \ 

Universal Gas Constant 

ytnole ■ K 

Reynolds number 

turbulent Reynolds number based upon integral length scale 

time dependent radius of laminar flame ball 

time dependent rate of increase of radius of laminar flame ball 

time dependent radius of turbulent flame ball 

time dependent rate of increase of radius of turbulent flame ball 

absolute flame speed 

consumption flame speed 

displacement flame speed 

displacement speed on the unbumed gas side 
initial energy deposited by spark 
laminar burning velocity 

planar unstretched laminar burning velocity 

laminar absolute speed 

laminar consumption speed 

laminar displacement speed 

turbulent burning velocity 

turbulent absolute speed 

turbulent consumption speed 

turbulent displacement speed 
time 

chemical time 
duration of spark 
turbulence time 

Kolmogorov time 
temperature 

Favre averaged temperature 
initial temperature 

adiabatic flame temperature 
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"A 

Uj 

u' 

U 

u 

Y a 

Y F 

Y F St 

Y 

p 

W 

x(t) 

Greek 


activation temperature 
velocity 

root mean square (rms) fluctuating turbulent velocity 

mean gas velocity 

mean velocity 

species mass fraction 

fuel mass fraction 

stoichiometric fuel mass fraction 

product mass fraction 
molecular weight 
position of particle j at time t 

Definition 


a 

a R 

Pr 

8 

L 

8 

T 

8(¥-<D 

Ar 

At 

AR 

£ 

V 

V 

sp 

A 

M 

v 

n 

r 

i//') 

p 

P h 

Pu 

<t> 

0 

® a 
OW 


specie 

reduced heat release factor 
reduced activation temperature 
laminar flame thickness 
turbulent flame brush thickness 

dirac function 

cell size in radial direction 

time step 

distance traveled by fluid particle in on time step At 
rate of turbulent energy dissipation 
Kolmogorov length scale 
spark efficiency 

Taylor microscale of turbulence 
dynamic viscosity 
kinematic viscosity 
mixing frequency 

reactive scalar vector for particle j 
density 

burned gas density 
unbumed gas density 
fuel equivalence ratio 
reduced temperature 

mass rate of addition per unit mass of species a due to reaction 

local fuel reaction rate -rate at which the reactants are consumed by the chemical 

reaction 
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Subscripts Definition 


b 

F 

i 

L 

sp 

st 

t 

0 


burned 

fuel 

particle i, coordinate direction 

Laminar 

spark 

stoichiometric 

turbulent 

initial 


Superscripts Definition 


a 

b 

c 

d 

j 

6 


66 


absolute 
modified 
consumption 
displacement 
particle counting 

root mean square fluctuating parameter, first derivative with respect to time 
flux, per unit area 


1. Introduction 

This paper is concerned with the numerical computation of early flame propagation for a spark- 
ignited, premixed methane/air gas in a moderate to intense turbulent environment. Real 
combustion processes occurring in such devices as internal combustion engines, rockets, natural 
gas burners, industrial power plants and gas-turbine engines nearly always occur within a flow 
field of moderate to intense turbulence. Several developments and modifications to the original 
algorithm have been implemented including a revised chemical reaction scheme and the 
evaluation and calculation of various turbulent flame properties. Several complex phenomena 
and topics are considered including spark ignition, turbulence, reacting flows and statistical 
methods. Attention is focused on the early stage of the flame front development as this stage is 
greatly influenced by turbulence and, thus, significantly influences overall engine performance. 
However, the aerodynamics of the flow, both convection and turbulence, is not the only 
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influence on the flame. As a result, vital influences such as the composition of the mixture and 
the characteristics of the ignition system such as energy delivered, spark plug geometry, and 
spark duration are also taken into account. 

Although numerous experimental and numerical studies on turbulent combustion have been 
developed, a common approach has not emerged due to the complexity of the physical and 
chemical phenomena occurring in turbulent flames. Recent reports [19, 20] have established the 
fact that a great deal of discrepancy still exists concerning how turbulent flame speeds are 
defined. The present work attempts to clarity these discrepancies by incorporating the absolute, 
displacement and consumption flame speeds as defined by Poinsot and Veynante [20]. 
Computing and accounting for the flame thickness and stretch and their effects upon the value of 
the flame speed will also be evaluated. Turbulent cases focus on the effect of varying turbulence 
intensity upon such characteristics as flame position, turbulent flame speed, and turbulent flame 
regimes. Attention is focused upon the distributed reaction zone where moderate to intense 
turbulence conditions exist. Traditionally, the theoretical foundation of the distributed reaction 
zone has been evaluated differently than the wrinkled flame regime and the corrugated flame 
regime. Recently within the turbulent combustion community, particularly in experimental 
research [7], there have been conclusions drawn that the distributed reaction zone is not as 
vulnerable an environment as traditionally believed and that a viable flame front can be sustained 
without the imminent threat of flame quenching. 

Computations involving turbulent conditions focus on the effect of parameters such as the 
ratio of the root mean square (rms) turbulent velocity to the laminar burning velocity, u '/£ , and 

the turbulent Reynolds number, Re v on various flame properties. Parameters such as the 
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ignition energy required and the quenching distance are verified. The described method 
produces results similar, if not identical, to those produced experimentally and numerically. 

2, Formulation of the problem 

In the present study chemical parameters have been set to simulate various premixed methane/air 
combustible mixtures at nonnal temperature and pressure (298.15 K and 1 atm). Lean 
methane/air flames are investigated with fuel equivalence ratios ranging from 0.6 to 0.9. The 
premixed gaseous mixture is initially quiescent in a non-decaying turbulent environment with the 
turbulence intensity varying from 0.5 m/s to 2 m/s. The mixture is assumed to be a perfect gas 
and one-step Arrhenius chemistry is utilized. A one-dimensional, spherically symmetric 
geometric model is considered and is shown schematically in Figure 1. The calculations 
simulate ignition of the premixed gaseous mixture by a spark ignition device. If the appropriate 
thermal requirements are met a spherical flame propagates radially outward from the point of 
ignition. The present analysis focuses on the flame front and flow conditions within 
milliseconds after successful ignition. During this time the combustion chamber pressure does 
not significantly change [16]. Asa result, a constant pressure environment at 1 atm is assumed. 

All results and calculations are made with respect to the center of the initial spark. The 
random motion of the center of the flame ball is not taken into account. Thus, results, such as the 
change in flame radial position, are calculated relative to the center of the initial spark. 

Finite difference solutions of the Navier-Stokes equations are often not practical, mainly 
because of the large dimensions of the PDF's. Thus, in the present study the approach chosen for 
obtaining a solution is an equation for the joint PDF of chemical and thermodynamic properties 
in turbulent reactive flows. The transport equation for the joint PDF is modeled and solved by a 
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Monte-Carlo method that simulates the representative terms in the PDF transport equation by the 
method of fractional steps. A Gaussian distribution for the normalized fuel concentration and 
fluctuating velocity is presumed. In the Monte-Carlo particle tracking method the fluid particles 
advance in space and are taken through representative collisions and chemical reactions. Each 
particle’s position coordinates, velocity components, and internal states are stored and modified 
with time. 

2. 1 Probability Density Function transport equation 

Starting from the conservation equations, with the assumption of low mach number and high 
Reynolds number, the properties at any point within the flame can be determined from a PDF 
transport equation [21]. Pope [22] proposed the use of a transport equation for the joint PDF of 
velocity, viscous dissipation, and reactive scalars. However, this equation does not include 
scalar gradients and therefore contains no information about the mixing time scale. Dopazo [9] 
therefore advocated the use of a transport equation for the joint statistics of velocity, velocity 
gradient, reactive scalars, and their gradients. In this equation, in addition to convection and 
chemical reaction being closed, so also is the tenn that describes the straining and rotation of 
scalar gradients, a mechanism that is believed to be essential in turbulent reacting flows. The 
closure problem is, however, shifted to the mixing of scalar gradients. This formulation has not 
yet been applied to flows with combustion [17]. 

In the current analysis the transport equation for the joint PDF of velocity and reactive 
scalars is considered [22]. The reactive scalars, such as the mass fraction of reacting species (the 
progress variable), can be represented by the vector i// The PDF transport equation provides a 

description of the transport of the PDF /?(\|/) in position and composition (species) spaces. The 
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main interest of the PDF transport equation is that the chemical reaction tenn depends only on 
chemical variables and does not require any modeling. Thus, the effect of chemical reaction on 
transport in composition space appears in closed form. 

The transport equation for the density-weighted joint PDF is: 


~ - Dp(.W ) + dpu'pXv) 

Dt dx t 

MeanFlowConvection TurbulentTransport 


dpp{y/)(b a {y/) d d r a^ a 

dy a dYa ^ dx i dx i 

Chemical Re action MolecularMixing 


( 1 ) 


where P '( M 7 ) is written for Si'P-cf)). The derivation and in-depth explanation of this equation is 
provided by Pope [21]. The terms on the left-hand side of equation (1) represent convection by 

the mean flow and by turbulence; they account for the transport of p(i/ z) in physical (i.e. flow 

field) space. The terms on the right hand side account for the transport of p(ip) in composition 
space due to chemical reaction and mixing. For chemically reacting flows it is of particular 
interest that the chemical source terms can be treated exactly for arbitrarily complex chemical 
kinetics. It has often been argued that in this respect the transported PDF formulation has a 
considerable advantage compared to other formulations. For chemically reacting flows the last 
tenn in (1), which represents the transport in reactive scalar (composition) space by molecular 
fluxes, is the most difficult to model. This term represents molecular or micro-scale mixing. 

Although the velocity-composition joint PDF is a function of many independent variables, 
its evolution equation can be solved by a Monte-Carlo method. Virtually all numerical 
implementations of PDF methods for turbulent reacting flows utilize Monte-Carlo simulation 
techniques [21, 17] in which the joint PDF is represented indirectly by a large number of 
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particles. The Monte-Carlo method tracks the individual particles through each of the 
phenomena represented in equation (1) by the method of fractional steps. 

2.2 The Monte-Carlo particle tracking method 

The Monte-Carlo solution strategy utilized is based on the concepts presented by Pope [22]. 
Monte-Carlo methods employ a large number of particles that can be considered as different 
realizations of the turbulent reactive flow problem under investigation. In the Lagrangian 
algorithm [21] the particles are not bound to grid nodes. Instead, each particle has its own 
position and moves through the computational domain with its own instantaneous velocity. 
Thus, the state of the particle is described by its position and velocity, and by the values of the 
reactive scalar that its represents as a function of time. The simplest method to estimate local 
means is to compute cell averages. Computationally, the radial direction is composed of a 
prescribed number of "cells" each with a specified radial length, Ar. Each cell is initially 
assigned a certain number of fluid particles with each particle having a prescribed velocity, fuel 
concentration, and temperature. For the results reported herein, each cell initially contains 1000 
particles. As time marches forward the number of particles within each cell changes as the flame 
propagates radially outward. 

In the Lagrangian simulation, typically the method of fractional steps [21, 17] is used. This 
method is based on the observation that the various terms describing the time evolution of the 
PDF in equation (1) are additive. Therefore the processes in physical, velocity, and reactive 
scalar space may be treated sequentially rather than simultaneously. For one-dimensional 
calculations the Lagrangian motion of particle j is given by 
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or 


dx 0) 

dt 


u^\x). 


( 2 ) 


x 0) (t + At) = x 0> (t) + u 0) A t 


(3) 


where the velocity u (i> describes the motion of particle j. According to the theory of stochastic 

differential equations the evolution of the spatial distribution of particles is represented by the 
first tenn on the left hand side of the transport equation for the density weighted joint PDF 
(equation 1). Similarly, the change of the composition value y/i of the reactive scalars of the 
Lagrangian particle j, is given by 


P 


n d ¥ 


<i) 


(j) 


dt 


= co 0) . 


(4) 


or 

c 

\\if(t + At) = \\ff(t) + ^yAt. (5) 

p J 

Equation (5) represents the first term on the right hand side of the PDF transport equation. This 
illustrates how the solution of only time-dependent Lagrangian equations for particles simulates 
the solution of the transported PDF equation. 

The main drawback of Monte-Carlo methods is that they suffer from a statistical error 
that decreases only slowly with the number of particles N pc per cell [17]. The error is 
proportional to 1/ A pr . For an acceptable numerical accuracy far more than a hundred 
particles must be present in each cell. For industrial CFD problems, which require large 
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numerical grids consisting of typically several hundred thousand grid cells, this leads to very 
large numbers of particles. Xu and Pope [26] have quantified the different errors that occur in 
Monte-Carlo methods for turbulent reactive flows. Monte-Carlo techniques are known to be 
sensitive to the number of cells, the number of particles, the mixing model, and initial conditions 
(i.e. spark shape and duration). 

3. Numerical procedure 

The numerical procedure is a fluid particle tracking method that approximates equation (1). The 
convective transport of the particle by the mean flow and by turbulence is accounted for in space 
and in time by the tenns on the left-hand side of the transport equation. The chemical and 
thennal composition of each particle changes due to molecular mixing and chemical reaction as 
represented by the terms on the right-hand side of the transport equation. 

3. 1 Input parameters 

Due to the large number of parameters involved in a spark ignited, turbulent flame analysis 
adjustment of parameters can enhance numerical results. It is important to note that in the 
current modeling, no ad-hoc adjustment of any parameter has been made. All parameters have 
their correct values throughout the analysis and were obtained from experimental measurement, 
theory or from the results of a thermodynamic analysis. When comparison is made to the work 
of other researchers, all ignition, mixture and flow conditions reported by the researcher are 
implemented in the present work. 

The following input values are required: time step, At, cell width, Ar, turbulent kinetic 
energy, k, dissipation of turbulent kinetic energy, e, normalized mean fuel concentration, c p , 
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root-mean-square fluctuation of fuel concentration, c'p , chemical time, t c h em , coefficient of 
ignition, S ign , pressure, P, spark duration, t sp , spark gap, L sp , activation temperature, T, u adiabatic 
flame temperature, Tf, initial reactant temperature, T 0 , pre-exponential factor, B, and the laminar 
flame speed, Sl. 

A chemical equilibrium calculation with dissociation using CEC Gordon and McBride code 
[11] is utilized to determine the adiabatic flame temperature. The adiabatic flame temperature is 
required in the calculation of the reduced temperature, 6, and in the Arrhenius reaction rate 
model. Table 1 provides values of the adiabatic flame temperature predicted by the algorithm. 


Table 1 Adiabatic flame temperatures for Methane/Air as computed 


using CEC [11] for lean mixtures as a function of (f> 


Methane/ Air Fuel 

Adiabatic Flame 

Equivalence Ratio, (j) 

Temperature, T f ( K) 

0.6 

1554.76 

0.7 

1718.25 

0.8 

1871.33 

0.9 

2011.50 


3.2 Chemical time 

The chemical time is defined as t c h em = v/Sl where v is the kinematic viscosity. The chemical 
time is required in the calculation of the turbulent Damkohler number, Dar- Various values of v 
and Sl are reported in combustion literature. Kinematic viscosity values of 15.5 x 10' m /s, 15.7 
x 10' m /s and 16.0 x 10' rn/s for a methane/air mixture are reported by Giilder [12]. The 
laminar burning velocity is dependent upon the fuel equivalence ratio, (j). Values of Sl can vary 
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based upon the various formulations including numerical, experimental and asymptotic methods 
[19]. Therefore, the value of t c i,em can vary significantly depending upon the values of v and Sl 
utilized. Because the values of t c hem and Dar are very sensitive to the laminar burning velocity, 
Sl, reliable Sl data for lean premixed methane/air flames is important in the present study. This 
is especially critical for (f> < 0.7 where Sl becomes very small (~ 0.1 m/s) and the scatter in 
theoretically calculated Sl data can be as high as 200 percent [6], Throughout the present work 
the kinematic viscosity is assigned a value of 15.7 x 10' m /s as reported by Giilder [12] for a 
lean methane/air mixture. The asymptotically based values of Sl provided by Poinsot and 
Veynante [20] are implemented in the present study. Table 2 provides the calculated values of 
tchem and the values of Sl utilized in all Monte-Carlo simulations. 

Table 2 Chemical time and laminar burning velocity values utilized in the Monte-Carlo method. 


Methane/ Air Fuel 
Equivalence Ratio, (/) 

Laminar Burning Velocity 
Sl (m/s) 

Chemical Time 

tchem (/MS) 

0.6 

0.13 

0.929 

0.7 

0.21 

0.356 

0.8 

0.28 

0.200 

0.9 

0.36 

0.121 


3.3 Activation temperature 

Poinsot and Veynante [20] present the parameters used in flame speed computations for a 
methane/air flame modeled by a single step reaction. An activation energy of 83,600 J/moIe was 
utilized as it is in the present work. The activation temperature, T A = Ea/R u , utilized is 10,055 K. 
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3.4 Turbulence parameters 

The constants calculated in the algorithm include the integral length scale of turbulence, l p 
turbulence time, t p and the turbulent Reynolds number, Re p which are defined as follows: 


r 3 


l T = — =C lr 


K 


1.5 


t T = ~ 

1 U 


Re , 


u'L 


(6) 

( 7 ) 

( 8 ) 


The above definitions for turbulence time and turbulent Reynolds number are generally 
agreed upon. However, the coefficient C, of the integral length scale of turbulence assumes 

various values throughout literature [1, 18, 23] depending on the exact method of formulation. 
Isotropic turbulence is assumed in the present work. Therefore, k = (3/2 )u' 2 or u' = ^J(2/3)k 
and the integral length scale of turbulence and the turbulence time are defined as: 



f 2 T ' 5 

Jv 



( 9 ) 


t 


T 



( 10 ) 


15 



The turbulence is non-decaying. Thus, the turbulence kinetic energy and the dissipation of 
turbulent kinetic energy are constant throughout a computation. This assumption does not 
introduce large errors since the computational time is typically 5 ms and the turbulence times are 
shown in Table 3 range from 0.121 ms to 0.929 ms. 

3.5 Gaussian distributions 

The initial fuel concentration and fluctuating velocity are prescribed a Gaussian distribution. 
After the initial Gaussian distribution is calculated for the fuel concentration, the subsequent 
concentration distribution proceeds randomly in time changing only due to mixing and reaction. 
The Gaussian distribution for the fluctuating velocities is updated at every interval of the 
turbulence time. 

The fluctuating velocities have a spatial correlation that is a function of the turbulent 
length scale. The correlation follows a simple zero order Frankiel decreasing exponential 
calculation. The random deviates of the Gaussian distribution are calculated using the Box 
Muller method. The normal deviates have a zero mean and unit variance. Initially each particle 
is assigned a fluctuating velocity and a fuel concentration value. Both the mean and root mean 
square fluctuations of velocity and fuel concentration are specified. Initially, since the premixed 
gas is at rest, the mean velocity is zero. 

3.6 Particle diffusion 

During each time step the particles in each cell move a certain distance based upon the 
instantaneous velocity and time step, At. The instantaneous velocity, u, is the sum of the mean 
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velocity, U, and the fluctuating velocity, u\ Thus, during each time step, each cell either loses or 
gains particles. For example, in the one-dimensional analysis the distance, AR, traveled by a 
fluid particle is defined as follows: 


AR = (U+u')At (11) 

Initially, the mixture is stationary and U(Q) = 0. These calculations provide the simulation of the 
propagation of the flame. As the flame propagates outward from the point of ignition, the density 
of the particles behind the wave decreases. Simulation of wave propagation is imperative since 
successful ignition of most combustion systems depends not only on the ignition of the fuel- 
oxidizer mixture, but also on the ability of the spark ignited flame kernel to propagate from the 
point of ignition. In many combustion systems it is in fact flame propagation rather than 
ignitability which detennines the flammability limits. 

3 . 7 Micro-scale mixing 

These calculations follow the model of Curl [8]. More elaborate mixing models exist; however, 
at the expense of increased computational requirements and complexity. Curl's model is a 
coalescence and redispersion model of the concentration and temperature. A random number 
generator randomly selects colliding particles. For example, when two particles of concentration 
Cj and C 1 mix one particle is formed with a concentration C = C 1 + C 2 . Upon redispersion there 

are two particles each with a concentration C/ = C 2 '= — (c, + C 2 ) . During small scale mixing not 

all of the particles will collide. The number of particles that mix is determined by the mixing 
frequency. The mixing frequency, !?, is the percentage of particles mixed per each time step. 
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The following empirically derived mixing frequency utilized by Borghi and Mantel [4] is 
employed in the present analysis: 


n= 0.625 



( 12 ) 


3.8 Chemical reaction 

The chemical reaction is represented by a single-step Arrhenius chemistry mechanism. 


R (reactants) P (products) 


(13) 


and the reaction rate co R is expressed as, 


cb R = B p c F exp 


(-C 



(14) 


where B is the pre-exponential factor. This expression is more conveniently cast in the fonn 
following [25,18,20]: 


cb R = B p c F exp 


' ziA 

y \-a R {\-0) y 


(15) 


©is the reduced temperature, 
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T-T 

<9=—^, (16) 

1 f 1 o 

where T o is the unbumed fresh gas temperature and T f is the adiabatic flame temperature for 
unity Lewis number. The activation temperature is T A , and the coefficients B , c/r, and J3r are, 

respectively, the reduced pre-exponential factor, the reduced heat release factor and the reduced 
activation temperature, defined as follows: 


B = B exp 


-Pr 

OCr 


( 17 ) 


cc R 



(18) 


Pr a R 



(19) 


The two parameters cxr and J3r measure the heat released by the flame and the activation 
temperature respectively. The values of (Xr and ()r used in the present work for T t = 1 0.05 5 /f are 
presented in Table 3. 
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Table 3 Values of ccr and J3r used in the Reaction Rate Model 


Methane/ Air 

Adiabatic Flame 



Fuel Equivalence Ratio, (/) 

Temperature, Tf(K) 

OCR 

Pr 

0.6 

1554.76 

0.8083 

5.23 

0.7 

1718.25 

0.8266 

4.84 

0.8 

1871.33 

0.8407 

4.52 

0.9 

2011.50 

0.8519 

4.26 


3.9 Spark ignition 

The spark energy decreases linearly in time and decreases parabolically with distance from the 
point of ignition [5, 6]. The spark does not produce any more energy after the specified ignition 
time, t , nor does it deposit energy at a distance exceeding the spark gap length, L . The total 
energy deposited by the spark in the spherically symmetric environment is [6]: 

7T 2 3 

^ign _ ^P^ign *sp ^ sp ‘ ( 20 ) 


The present ignition model contains a spark efficiency, p sp , to account for losses in actual spark 
ignition systems. In actual spark ignition systems much of the electrical energy is expended in 
radiative losses, shock wave formation, and convective and conductive heat losses to the 
electrodes and flanges [10]. Flanges connected to the electrodes are often used in actual spark 
ignition systems. Zeldovich [27] reports the spark efficiency for mixtures varies from 2-16%. 
The spark efficiency is the ratio of external input energy required in order to heat a spherical 
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volume of gas to its adiabatic flame temperature to the total energy in a capacitance spark. 
Therefore, the energy required to ignite a mixture is equal to 2-16% of the electrical energy in 
the spark due to the losses described by Glassman [10]. In the present computer model a spark 
efficiency of 10% is employed in order to computationally account for and model actual losses 
associated with real world spark electrodes. 

3.10 Average cell temperature and density 

During each time step the Favre averaged temperature of each cell is calculated as follows: 

n 

Y.PiTi 

T = i=l ( 21 ) 

n 

YjPi 

1=1 


where n is the number of particles in the cell during the current time step. Upon substitution of 
the ideal gas law in equation (21), the average temperature can be expressed as follows: 


T = 



I 


i=l 


P 


RTi 


or 


T = 


n 


" 1 

If 

i=\ 1 


(22) 


(23) 
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The Favre averaged density of each cell is calculated using the ideal gas law at constant 


pressure, P, and the computed Favre average temperature, T , of each cell: 


P = 


P 

R-u y 
w 


(24) 


where W is the mean molecular weight of the mixture given by: 


W 


N Y 

Z Y a 

W r ' 


a= 1 


a 


(25) 


Initially, for the unburned fresh gas methane/air mixture, W = 29 kg/kgmol. Therefore, for the 
unbumed gas mixture, equation (24) can be written as: 


P 

P = = 

2S7(J/kg-K)T 


(26) 


The reduced, non-dimensional temperature of each cell is calculated using equation (16). 


3.11 Average cell velocity 

During each time step the discretized equation of continuity is solved for the average velocity of 
each cell: 
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1— A (pUr 1 ) ^ ( 27) 

A? ? .2 Ar 

Initially, the gas mean velocity is zero or U (t = 0) = 0. pU is the ensemble or time averaged, 
density weighted velocity at a fixed radial location. 

3.12 Flame front position 

The numerically simulated results are analyzed similarly to other numerical and experimental 
measurements [2, 23] that base the flame position on the point where the normalized fuel 
concentration, or temperature, is the average of the hot and cold boundaries. 

3.13 Post-processing 

The Gordon and McBride CEC code [10] is also utilized to verily the burned to unbumed gas 
density ratio used in the calculation of the turbulent burning velocity. Aung et al. [2] similarly 
utilized the CEC algorithm to determine the density ratio. It should be noted, however, that this 
approach represents a convention that ignores preferential diffusion effects. These effects 
modify the local mixture ratio and thennal energy transport for stretched flames, and thus the 
local density ratio of the flames. However, as noted by Aung et al. [2] this convention is 
convenient because a single density ratio is used to relate flame speeds (i.e. absolute speeds) and 
burning velocities or displacement speeds for all levels of flame stretch. This avoids the current 
uncertainties about the effects of stretch on the jump conditions across flames for particular 
conditions in the unburned gas mixture. Numerical simulations provide a means of estimating 
density ratio variations but there are still uncertainties concerning the reliability of the 
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calculations. Furthermore, Aung et al. [2] estimate that variations of flame properties due to 
variations of density ratios with stretch do not exceed 11%. This variation is significantly less 
for fuel lean conditions where Markstein numbers, and thus the effects of stretch, are small. This 
contrasts fuel-rich conditions where Markstein numbers, the effects of stretch, and thus, the 
variation of properties due to density ratio variations are large. 


Table 4 Burned to Unburned gas density ratio computed using 
CEC code [10] as a function of $ 


Methane/ Air Fuel 

Density Ratio, 

Equivalence Ratio, (j) 

pi/ pll 

0.6 

0.1859 

0.7 

0.1675 

0.8 

0.1534 

0.9 

0.1420 


4. Preliminary results 

4.1 Monte-Carlo sensitivity 

As previously mentioned, Monte-Carlo methods are sensitive to input parameters such as the 
spark ignition characteristics and the number of particles initially assigned per cell. A 20 percent 
difference in calculated flame property values can result when 100 particles per cell are used 
instead of 1000 particles per cell. A plot of the flame front radial position versus time from 
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spark onset is provided in Figure 2. The results are based upon the initial assignment of 100 
particles per cell. Figure 2 provides the results of 10 identical turbulent flame simulations using 
the present Monte-Carlo method. The average curve of these ten simulations is also presented. 
Figure 3 is based on simulations identical to Figure 2 except that 1000 particles were initially 
assigned to each cell. It is obvious that utilizing 1000 particles per cell generates less scatter in 
the results. The results of individual executions using 1000 particles per cell are much more 
reliable than individual executions generated using 100 particles per cell. Also, the average 
flame radial position generated utilizing 1000 particles travels faster than the flame generated 
using 100 particles per cell. At approximately 8 ms the flame position in Figure 2 is 
approximately located at 4.3 mm. In Figure 3, at 8 ms, the flame radial position is located at 
approximately 5.0 mm. Thus, there exists approximately a 16 percent difference in flame radial 
position between 100 particles per cell versus 1000 particles per cell. In the present work 1000 
particles are initially assigned per cell. All results provided are based on the averages of 10 
individual simulations in order to reduce error. 

4.2 Quenching distance 

The quenching distance is the critical distance between spark electrodes that will allow an 
inflamed spark kernel to grow unaided. Experiments conducted by Lewis and Von Elbe [15] 
utilized flanges connected to the tips of the electrodes. These flanges were made from an electric 
nonconductor such as glass so that the sparks remained centered between the electrodes. 
Calculations for (f> = 0.9 (methane/air), with u'= 0.0 m/s and c' F = 0.0, using the current 

numerical model show the quenching distance to be 2 mm. This matches exactly the 
experimentally determined quenching distance of Lewis and Von Elbe [15] for the same mixture. 
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The quenching distance is determined by varying the spark electrode gap distance, L , while 
keeping all other conditions and variables constant until the point of ignition is reached. For 
values of L sp up to 1.9 mm the flame does not ignite. As L is increases from 1.9 mm to 2.0 mm 
the flame ignites and propagates outward from the center of the initial spark. 

4.3 Minimum ignition energy 

The minimum ignition energy is the least amount of energy required from the spark to create a 
spark kernel at the quenching distance. The current model predicts a required minimum ignition 
energy, S jan mjn , of 0.25 to 0.3 mJ for L sp = 2.0 nun and <j>= 0.9 at 1 atm. This value was 
determined in the same manner as was the quenching distance. The initial energy deposited by 
the spark, S ign , was increased, while all other conditions and variables remained constant, until 
ignition and flame propagation was observed. This compares very well to the experimental data 
of Lewis and Von Elbe [15] who report a required minimum ignition energy of 0.28 mJ to ignite 
an equivalent mixture of methane/air for L sp = 2.0 nun and (f> = 0.9 at 1 atm. Blanc et al. [3] also 
report a minimum ignition energy of 0.28 mJ for the same mixture. 

5. Turbulent flame results 

The purpose of the present results is to demonstrate that the current model is a viable predictor of 
the effect of the mixture composition, the aerodynamics of the flow and the characteristics of the 
ignition system upon early flame propagation in a moderate to intense turbulence environment. 

The turbulent flame speed and burning velocity is not only a function of the fuel-oxidizer 
mixture ratio but also of the flow characteristics and the ignition source configuration. 
Therefore, there can exist great difficulty in correlating both experimental and numerical data of 
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various investigators. When comparing the present results to various experimental and 
numerical results of other researchers, every attempt has been made to simulate the conditions of 
the referenced work. 

The exact fonnulation of all Monte-Carlo results with respect to the turbulent absolute, 
displacement and consumption speed will be established. When comparison is made to the 
flame speed results of other researchers and no definitive clarification was provided in regard to 
a kinematic or reaction based fonnulation, the turbulent flame speed is referred to as the 
turbulent burning velocity. However, displacement based fonnulations are most commonly 
reported in the literature even though not directly specified. When the flame thickness is 
neglected the displacement and consumption speeds are approximately identical. Evaluation of 
the turbulent flame regime criteria, the effect of turbulence intensity and predictions of various 
turbulent flame speeds calculated using the present Monte-Carlo code are presented and 
compared to the values of other researchers. Attention is focused on flames in the distributed 
reaction zone regime. 

5 . 1 Turbulent Flame Properties 

Flames with varying degrees of turbulence intensity are examined. The properties associated 
with these flames are shown in table 5. In each case the unburned gas temperature is 298 K, the 
activation temperature is 10,055 K, the pressure remains constant at 1 atm and the kinematic 

viscosity, v, is 15.7 x 10' m /s [12]. Note the following notation utilized in Table 5: laminar 

displacement speed, S^, laminar flame thickness, Sl, modified laminar flame thickness, S b L , 

Kolmogorov time, t n , Karlovitz number, Ka, Karlovitz flame stretch factor, Ka ' , Kolmogorov 
length scale, 77 , Taylor microscale of the turbulence, A. 
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5 .2 Turbulent Flame Regimes 

Classification of turbulent premixed flames continues to evolve after various regimes were 
originally defined over 60 years ago. The classification presently utilized follows the standard 
prescribed by the Borghi-Barriere diagram [7]. In this diagram the Ka = 1 boundary, commonly 
referred to as the Klimov- Williams criterion, separates the regimes of wrinkled laminar flames 
from flames with thicker distributed reaction zones. Currently, among combustion researchers, 
there is much interest in verifying experimentally the changes in flame structures as predicted by 
the Klimov-Williams criterion [7]. This is due to the fact that many practical systems span the 
criterion boundary and theoretical treatments of the two flame regimes can be quite different. 
For values of Ka > 1 the flame is defined as being in the distributed reaction zone regime. 
Traditionally the distributed reaction zone regime has been viewed as incapable of sustaining a 
flame due to the fluctuations in temperature, concentration and velocity. Recently however, 
numerical and theoretical data seems to suggest that flamelets are more resilient to penetration by 
small eddies than is prescribed by the Klimov-Williams criterion. 

Monte-Carlo results for various flames seems to support the latest views that the 
distributed reaction zone regime of Ka > 1 should be interpreted as the regime where the 
probability for small intense turbulence to penetrate the flame sheet is finite [7]. Based upon the 
Borghi-Barriere diagram most of the turbulent flames currently simulated are specified as 
occurring in the distributed reaction zone regime. Table 6 provides specific values of regime 
criteria as a function of the laminar flame thickness, 8l, and a modified laminar flame thickness, 
S h L . The modified laminar flame thickness is a function of the adiabatic and initial flame 
temperature and is a more precise estimate of the laminar flame thickness [20], 
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Table 6 Flame Regime criteria 


Flame 

r 

U 

"Si 

Sl 

s\ 

Ka 

Flame 

Regime 

FI 

2.5 

25 

3.5 

0.791 

Corrugated 

F2 

2.5 

52 

7.3 

1.40 

Distributed 

F3 

2.5 

48 

6.7 

1.53 

Distributed 

F4 

2.5 

50 

6.9 

1.57 

Distributed 

F5 

2.5 

52 

7.2 

1.64 

Distributed 

F6 

2.5 

54 

7.4 

1.69 

Distributed 

F7 

2.5 

55 

7.6 

1.75 

Distributed 

F8 

2.5 

57 

7.9 

1.80 

Distributed 

F13 

2.5 

69 

9.5 

2.17 

Distributed 


Depending upon which value of the laminar flame thickness is used, 8l or 8 h L , a large difference 
exists in the values of the ratio of the integral length scale to the laminar flame thickness. This 
difference is readily apparent in Figure 3 in which flames F1-F8 and FI 3 are assigned to their 
respective flame regimes based upon the criteria presented in Table 6. 

5.3 Concentration, Density and Velocity Profiles 

Time evolution of the nonnalized fuel concentration, density, and mean velocity of a typical case 
run represented by the parameters associated with flame FI are shown in Figures 5, 6, 7, and 8. 
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Figure 5 displays the time history of the normalized fuel concentration versus the flame radial 
distance. It is readily observed that the flame is propagating radially outward. The value of the 
fuel concentration drops to zero indicating complete combustion at the center. It is from the 
nonnalized concentration profiles that the flame radial position and absolute flame speed are 
calculated. The radial position of a flame front is taken where the nonnalized concentration is 
the equivalent of half of the initial normalized fuel concentration. This is a common convention 
used by various combustion researchers [23]. 

The gas density plotted against the radial distance for various times from spark onset for 
flame FI is shown in Figure 6. The ratio of the burned gas density to the unburned gas density, 
pjp u , is approximately 0.1. The ratio predicted by the Gordon McBride CEC code [11] for ^ = 
0.8 is 0.1534. This value is common for most practical gaseous hydrocarbon flames in which the 
burned to unburned gas density ratio lies in the range 0.1 to 0.2 [1]. 

The mean gas velocity versus radial distance is shown in Figures 7 and 8. The large initial 
increase in the mean velocity obvious at 0.4 ms is a result of the initially rapid volume expansion 
created by the spark ignition energy. The spark ignition energy weakens due to an imposed 
linear decay in time and a parabolic decay in distance. The spark is ignited a duration of 1.5 ms. 
It can be observed in Figure 7 that between 0.4 and 1 .6 ms the mean velocity decreases with the 
lowest maximum values occurring at 1.6 ms. After the 1.5 ms spark duration the velocity again 
increases. 

Figure 7, in comparison to Figure 5 and Figure 6, shows that the highest radial gas 
velocities occur within the flame front. For example, it can be seen from Figure 6 that the 
maximum gas velocity at 2.8 ms occurs at approximately 3 mm. By comparison with Figure 5 
and Figure 6 it is observed that this position lies between the burned and unburned regions of the 
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flame front. Thus, the maximum mean gas velocity occurs within the turbulent flame brush. It 
has been well established that in all flames there is a large increase in velocity as the gases enter 
the burned gas state [10]. Therefore, heat release itself can play a role in inducing turbulence. 
However, the effect of heat release upon turbulence properties is not considered in the present 
work as turbulence is assumed to be non-decaying. 

5 . 4 Turbulent Flame Absolute Speed 

The turbulent flame absolute speed, Sj , is also referred to as the flame propagation speed or, for 
turbulent spherical flames, the rate of increase of the turbulent flame ball radius R' T ( t). Similar 
to the laminar cases, this speed can be obtained by calculating the slope of the time evolution of 
the flame radial position. Pope and Cheng [23] previously developed a Monte-Carlo numerical 
solution of the modeled equation for the joint PDF of the velocities and the reaction progress 
variable. Pope and Cheng [23] also present the turbulent flow experimental results generated by 
Hainsworth [14]. The experimental investigations generated turbulence by passing a premixed 
methane/air mixture through a turbulence-generating perforated plate producing a uniform mean 
flow of homogenous, isotropic turbulence. The mixture was ignited downstream of the plate 
forming a spherical flame kernel that expanded and was distorted by the turbulence. A 
comparison of the results of the present model with those of Hainsworth [14] and Pope and 
Cheng [23] is presented in Figure 8 which shows the time evolution of the flame front radial 
position. The turbulence intensity, u ' , is 1.93 m/s and (f> = 0.8. Other properties of the Monte- 
Carlo simulated flame are provided in Table 5 for flame FIS. As seen in Figure 9 between 0 and 
approximately 2 seconds S T is approximately equal to the laminar displacement speed, S d L . 
Between 2 and 5 seconds S T increases above the laminar value. This observation is consistent 
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with previous observations of Anand and Pope [1] and Rutland et al. [24]. Excellent agreement 
exists between the present Monte-Carlo method, the experiments of Hainsworth [14] and the 
statistical calculations of Pope and Cheng [23]. The absolute flame speed calculated by the 
Monte-Carlo method is 3.8 m/s. 


5.5 Turbulent Flame Displacement Speed 

Similar to the laminar displacement speed the turbulent displacement speed, Sf , is defined as 



( PdL 
Pu 


r't (0 - 


( BdL 

Pu 


o U 


(28) 


where R' r (/) is the rate of increase of the turbulent flame ball radius, R 7 (t). R\ (t) is the turbulent 
flame propagation speed and is also referred to as the turbulent flame absolute speed, S ^. . 
Recall, the results of the Gordon McBride CEC code [1 1] report a density ratio, p b l p u , of 0.1534 
for (f> = 0.8. Therefore, the present model predicts the turbulent displacement speed to be 
0.58 m/s. 

5. 6 Turbulent Flame Consumption Speed 

When the flame thickness is assumed to be negligible the turbulent displacement speed and the 
turbulent consumption speed are assumed to be equal. In order to obtain a more precise estimate 
of the consumption flame speed the flame thickness must be evaluated. This can be done by 
providing better estimates of the burned gas mass. The mass must include not only the fully 
burned gases but also the burned gases within the flame front itself. Poinsot and Veynante [20] 
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provide an expression for the consumption flame speed as a function of r(t) which includes the 
effects of finite flame thickness. For a turbulent flame the consumption speed is defined as 
follows: 



1 + + 

2 r(t) p b 


Pu dt 


(29) 


As specified in Table 5 for flame F13, 8 l = 56.01 pm and d b L = 405.3 pm. In order to correctly 
represent the impact of the turbulent flame consumption speed, calculation of the flame thickness 
is required as a function of time. Particularly for turbulent flames, the flame thickness increases 
in time. The turbulent consumption speed evaluated using 8l and a radial position of 10 mm is 
0.59 m/s. This varies little from the Monte-Carlo predicted displacement speed of 0.58 m/s. The 
evaluation of consumption flame speeds at larger radial locations is usually neglected due to this 
fact. Flowever, using the more precise modified laminar flame thickness, 8 b L , the consumption 
flame speed increases significantly to a value of 0.67 m/s. Therefore, the evaluation of the 
consumption flame speed based upon S b should not be neglected as the consumption speed 
varies greatly from the displacement speed even at larger radial distances. 

5. 7 Effect of Turbulence Intensity 

The radial position of the flame as a function of time from spark onset is shown in Figure 10 for 
seven flames with increasing values of turbulence intensity. The flames represented are F2, F3, 
F4 , F5, F6 , F7 and F8 for (f> = 0.8. The turbulence intensity ranges from 1.3 m/s to 1.6 m/s. 
The chemical and spark ignition conditions are identical in each of these seven flames. The only 
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variation occurred in the turbulent flow properties by increasing the turbulence intensity. The 
dissipation of the turbulent kinetic energy increased with each increase in the turbulence 
intensity in order to maintain a constant turbulence time of 2 ms in each flame simulation. Other 
properties associated with these flames are specified in table 5. It is observed that within the 
range of specified turbulence intensity values the turbulent absolute flame speed, S “ , and the 
turbulent displacement flame speed, S* , increases with u'. Initially, for approximately 2 ms, 
there exists little variation in flame propagation among the seven flames. The flame does not 
display any sensitivity to the value of the turbulence intensity and travels at approximately the 
value of the laminar flame speed for (f> = 0.8. At times greater than 2 ms flames F2-F8 
experience increasing rates of flame propagation. Turbulence convects, wrinkles, and stretches 
the flame sheet, increasing the flame area and thus the flame speed. 

6 Conclusions 

The effect of turbulence on a propagating premixed flame is a subject of long standing practical 
importance and continuing research. Incorporation of all influences into an ignition model and 
the use of appropriate turbulence models are essential for correct treatment of flame kernel 
formation. The present model considered several vital influences such as spark ignition energy 
losses, flow field effects, mixture properties and thennodynamic conditions upon early flame 
propagation. Emphasis is placed on the determination of various turbulent flame speeds 
including the absolute, displacement and consumption flame speed. 

The model predicts a minimum ignition energy, S ign _ min , of approximately 0.3 mJ for <j) = 
0.9. This agrees well with the experimental values of 0.28 mJ cited by Lewis and Von Elbe [14] 
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and Glassman [9]. The model predicts a quenching distance of 2 mm that also matches the value 
observed by Lewis and Von Elbe [14]. 

For a turbulence intensity of u' = 1.93 m/s the predicted flame ball position matches very 
closely to the experimental values predicted by Hainsworth [13] and the numerical calculations 
of Pope and Cheng [22]. The present model demonstrates that during early ignition the flame 
travels at the laminar burning velocity for approximately 2 ms then increases to an asymptotic 
turbulent flame speed value. The asymptotic turbulent flame speed computed for (j) = 0.8 is 
0.58 m/s. 

The results verify an increase in flame propagation speed with increasing initial turbulence 
intensity. It is shown that even for varying turbulence intensities the duration of the initial phase 
of the flame that travels at S L occurs for a period of approximately 2 ms. The present model 
predicts that increasing the turbulence intensity increases S T and that the duration of the initial 
laminar portion of turbulent flames is affected very little by variation of the turbulence intensity. 

Physically, the current model predicts the mean gas velocity to be highest inside the flame 
brush. There is significant radial motion of gases inside and well ahead of the flame in the 
unbumed gas region. It is observed that the turbulent flame brush thickness increases with 
increasing time from ignition. 

Calculation of the turbulent consumption speed based upon a modified laminar flame 
thickness provides a value of 0.67 m/s. This varies from the corresponding displacement speed 
by 15.5 percent and emphasizes the importance of the consumption flame speed. The effects of 
preferential diffusion cause the unburned-to-burned gas density ratio to vary as a function of 
flame stretch. However, it is difficult to calculate the variation of the density ratio across the 
flame front. Therefore, in the present analysis, a single density ratio was calculated using the 
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CEC Gordon McBride code [11]. Also, for fuel-lean conditions the effects of stretch are not as 
great as for fuel-rich conditions and, presumably, the density ratio does not vary significantly. 
This approach can be easily converted to correct the calculated burning velocities or 
displacement speeds once density ratios can be generated from well-established flame structure 
models. The theoretical treatment of flames in the distributed reaction zone was not handled 
differently than flames within the wrinkled or corrugated flame regimes. The Monte-Carlo 
results for various flames supports the latest views that flames in the distributed reaction zone 
regime are not as vulnerable to flame quenching as commonly believed. The distributed reaction 
zone regime has traditionally been viewed as incapable of sustaining a flame due to the 
fluctuations in temperature, concentration and velocity. 
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Figure Captions: 

Figure 1 . — Schematic of Spark-Ignited turbulent flame ball and Kinematic Balance between the 
Flame Propagation Velocity, Flow Velocity and Turbulent Burning Velocity. 

Figure 2. — Ten Realizations and Average of the radial flame propagation as defined by flame 
position as a function of time. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, 
u’ = 0.7 m/s, Initially Gaussian Distribution for velocity with U = 0.0 m/s. Initially 100 particles 
distributed per cell. 

Figure 3. — Ten Realizations and Average of the radial flame propagation as defined by flame 
position as a function of time. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, 
u’ = 0.7 m/s, Initially Gaussian Distribution for velocity with U = 0.0 m/s. Initially 1000 
particles distributed per cell. 

Figure 4. — Borghi-Barriere Regime diagram for turbulent premixed flames. 

Figure 5. — Time evolution of the normalized fuel concentration as a function of radial position 
for flame FI. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, u’ = 0.7 m/s, 
Initially Gaussian Distribution for velocity with U = 0.0 m/s. Initially 1000 particles distributed 
per cell. 
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Figure 6. — Time evolution of the cell averaged gas density as a function of radial position for 
flame FI. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, u’ = 0.7 m/s, Initially 
Gaussian Distribution for velocity with U = 0.0 m/s. Initially 1000 particles distributed per cell. 

Figure 7. — Time evolution of the cell averaged mean gas velocity as a function of radial position 
for flame FI. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, u’ = 0.7 m/s, 
Initially Gaussian Distribution for velocity with U = 0.0 m/s. Initially 1000 particles distributed 
per cell. 

Figure 8. — Three Dimensional Time evolution of the cell averaged mean gas velocity as a 
function of radial position and time from spark onset for flame FI. Initial flame properties: 
Premixed Methane/Air mixture, (j) = 0.8, u’ = 0.7 m/s, Initially Gaussian Distribution for velocity 
with U = 0.0 m/s. Initially 1000 particles distributed per cell. 

Figure 9. — Comparison of the flame propagation as defined by flame position as a function of 
time for flame FI 3. Initial flame properties: Premixed Methane/Air mixture, (j) = 0.8, u’ = 1.93 
m/s, Initially Gaussian Distribution for velocity with U = 0.0 m/s. Initially 1000 particles 
distributed per cell. 

Figure 10. — Effect of turbulence intensity on the flame propagation as defined by flame position 
as a function of time for flames F2, F3, F4, F5, F6, F7, and F8. Initial flame properties: 
Premixed Methane/ Air mixture, c|) = 0.8, u’ = 1.3 to 1.6 m/s, Initially Gaussian Distribution for 
velocity with U = 0.0 m/s. Initially 1000 particles distributed per cell. 
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